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Abstract 



This purpose of this write-up is to share an idea for accurate compu- 
tation of Laplace eigenvalues on a broad class of smooth domains. We 
[ J | represent the eigenfunction it as a linear combination of eigenfunctions 

-h ■ corresponding to the common eigenvalue p : 

^— > 

u (r, 9) = 2_^ PnJn (pr) cos nd, 



We adjust the coefficients P„ and the parameter p so that the zero level 
set of u approximates the domain of interest. For some domains, such as 
ellipses of modest eccentricity, the coefficients P n decay exponentially and 
the proposed method can be used to compute eigenvalues with arbitrarily 



,—4- ' high accuracy. 

in ' 

£D ' 1 Introduction 

The celebrated level set method |2| is a numerical method for solving problems 
with moving interfaces and has been applied to eigenvalue problems |3]- The 
method proposed here also represents interfaces as level sets of functions. flow- 
s' | ever, our level set function is a linear combination of global functions, with the 
V^ ■ coefficients of the linear combination and possibly additional parameters acting 

as the degrees of freedom for deforming the level set. When only a small num- 
ber of functions is required for the accurate computation of some quantity, the 
proposed method can be viewed as essentially analytical. 

The Laplace eigenvalues (with zero or any other type of boundary conditions) 
are not available analytically for a 2 x 1 ellipse. (In fact, eigenvalues for only a 
handful of shapes have been discovered since Raylcigh's pioneering "Theory of 
Sound" Hj.) Furthermore, eigenvalues are difficult to compute with high accu- 
racy. The finite element method [5] with isoparametric quadratic elements would 
perhaps be the most common approach. Combined with Richardson extrapola- 
tion, it is an impressive 5-th order method. However, in the author's experience, 
it fails to deliver the lowest eigenvalue beyond the 13-th digit. Furthermore, the 



finite element method is ineffective in computing the high (enough) eigenvalues 
since the supporting mesh must resolve the high frequency oscillations of the 
corresponding cigcnfunctions. 

2 Description of the approach 

We propose a method that computes the lowest eigenvalue on a 2 x 1 ellipse 
to arbitrary accuracy. The method can be applied to other shapes and higher 
eigenvalues. The method has a number of desirable features, including utmost 
simplicity, exponential accuracy, ease of analysis and the ability to produce an 
exact solution to an approximate problem. 

Represent the domain Q with boundary S as the interior Q u of the zero level 
set S u of the function u (r, 9) given by the linear combination 

N 

u (r, 9) = y^ P n J n (pr) cos n6>, (1) 

n=Q 

where J n are Bessel functions. The coefficients P„ and the parameter p are used 
as degrees of freedom in fitting the domain Q u to fi. Since u (r, 9) given by (QJ 
satisfies 

- Au = p 2 u, (2) 

then, by definition, p 2 is a Laplace eigenvalue on Q u with zero boundary condi- 
tions and u (r, 9) is the corresponding eigenfunction. Note that while Q u is an 
approximation to CI, the function u (r, 9) is an exact eigenfunction on Q u . 

The objective function measuring the proximity between S and S u can be 
defined as 

f D 2 u (S)dS, (3) 

Js 

where D u is an appropriately defined distance function that captures the dis- 
crepancy between S and S u . The Mathcmatica code below takes a pragmatic 
approach and defines D (a) as 

D(a)=r (a) - r u (a) , (4) 

where r (a) is the polar representation of the ellipse and r u (a) is the polar rep- 
resentation of S u . The minimization can be carried out effectively in a number 
of ways, including available generic optimization routines. 

As mentioned above, one of the most appealing characteristics of the pro- 
posed method is its utmost simplicity. The Mathcmatica code that computes 
the lowest eigenvalue on the 2x1 ellipse is given below in its entirety. 

1. Goals = { AccuracyGoal -> 75, PrecisionGoal -> 75, WorkingPrecision 
-> 120}; 

2. A = 1/2; B = 1; 

3. NumberOf Terms = 30; 

4. Shape [thetaj = 1/Sqrt [Cos [theta] ~2/A~2 + Sin[theta] ~2/B~2] ; 



5. lsf[rho_, p_] [r_, thetaj := BesselJ[0, rho*r] + Sum[p[[n]] BesselJ[2 
n, rho*r] Cos [2 n theta] , {n, 1, Length[p]}]; 

6. r[rho_, p_] [thetaj := x /. FindRoot [lsf [rho , p] [x, theta] == 
0, {x, Shape [theta] } , Evaluate [Goals] ] ; 

7. QbjectiveFunction[rho_, p : { ?NumericQ}] := Sum[(r[rho, p] [theta] 

- Shape [theta]) ~2, {theta, 0, 2 Pi - 2 Pi/60, 2 Pi/60}]/ (2 Pi/60) // 
Sqrt; 

8. s = FindMinimum[ObjectiveFunction[rr , pp] , {{rr, BesselJZero [0, 

!]}> {PP> ConstantArray [0 , NumberOf Terms] }} , Evaluate [Goals] , Maxlterations 
-> 500, StepMonitor :> Print [rr, " ", ObjectiveFunction[rr , pp] // 
N[#, 5] &, " ", Date[]]]; 

Code notes: 

Line 1 specifies accuracy goals and the working number of digits. Naturally, 
higher working precision takes more time. On a circa 2007 desktop, the pre- 
sented code takes about 30 minutes to take the first step and about 10 minutes 
for each subsequent step. It takes several steps to reduce the error by a factor 
of 10. These figures are given for computing with 30 terms in series ([!]) and 
125 digits of accuracy. The method is substantially faster when fewer digits are 
required. For instance, 6 digits of accuracy can be obtained with 3 terms in 
series {T]) in a matter of seconds. 

Line 2 specifies the semiaxes of the ellipse. 

Line 3 specifies the number of terms in equation (TTJ) . Due to the symmetry 
of the ellipse with respect to the y-axis, we only use the even terms in the series 
(fTJ), raising the effective number of terms to 60. The variable NumberOfTerms 
is referenced in Line 8. 

Line 4 gives the ellipse in polar coordinates 

Line 5 defines the level set function (lsf) u (r, ff) . This line specifies that 
u also depends on the parameters p and P n . Since the zero level set remains 
unchanged when the coefficients P n are multiplied by a number, the expansion 
here assumes that Pa = 1 and gives u (r, 9) as 

N 

u (r, 6) = J {pr) + ^ P 2n J 2n (pr) cos 2nd. (5) 

n=l 

Line 6 defines the zero level set in polar coordinates by solving the equation 
u (r, 6>) = 0. 

Line 7 evaluates the objective function in equation ([3]). The integral is 
approximated by a finite sum with 60 terms. More terms are needed when 
greater accuracy is targeted. 

Line 8 finds the optimal configuration. The initial configuration is the unit 
circle (since all P n = for n > 1). 



3 Further thoughts 

3.1 Hadamard acceleration 

The level set function u (r, 9) is an exact eigenfunction on Q u . This fact can be 
used increase the accuracy of the eigenvalue estimate. The Hadamard formula 
P] gives the rate of change in the eigenvalue A in response to a deformation of 
the domain tt. Consider a smooth evolution S (t) parameterized by t. Then the 
rate of change dX/dt in the eigenvalue is given by 

§ = -[ C\V^\ 2 dS, (6) 

at JS(t) 

where C is the Hadamard velocity of 5* (t) and ip is the eigenvalue corresponding 
to A. 

Since p 2 is an exact eigenvalue on S u and u (r, 8) is the corresponding exact 
eigenfunction, the Hadamard formula can be effectively used by letting C equal 
the distance between S u and S along the normal direction to S u . 

3.2 Prescribing C 

Deforming the surface according to a prescribed C is likely a more effective 
way of finding the optimal u than using a generic optimizer. Suppose that D u 
is the signed normal displacement between S and S u . Then the differential 
equation C = —D u on the moving interface S u corresponds to S u approaching 
S. Therefore, we may vary the parameters P n and p in such a way that the 
resulting velocity of S u is as close as possible to — D u . 

Let us calculate the velocity of the interface C that results when the pa- 
rameters of the level set function u are varied in a prescribed way. In a general 
setting, suppose that the function u depends on the parameters P n and we are 
now including p among the P n . Refer u arbitrary curvilinear coordinates Z±. 
Suppose that the parameters P n vary according to P n (t) and the resulting zero 
level set has the equation Zi (t,S), where S are arbitrary coordinates on the 
surface S u . The functions Zi (t, S) satisfy the equation 

u(P n (t),Zi(t,S)) = (7) 

Differentiating with respect to t, we find 

dF ■ dF dZ, (t, S) 

where P n = dP n /dt and summations over n and i are implied. Denote the 
gradient OF/dZi by |VF| N (where Ni is the unit normal) and dZi/dt by Vf. 

dF ■ 

—P n + |VF| NiVi = 0, (9) 



Since £ N.V, = C, we find 



C=-Pn^\*F\-\ (10) 



Therefore, P n should be chosen so that the series (|10|) approximates C as closely 
as possible - for example, in the least squares sense. 

3.3 Enriching the family of functions 

The presented calculation is based on a very limited set of functions J n (pr) e'" s , 
all centered at the same origin. Naturally, any Laplace eigenfunction (corre- 
sponding to the same eigenvalue) can be added to the mix. For the problem at 
hand, it is beneficial to consider functions of the form J n (pr) e , shifted other 
poles. Different geometries may utilize other functions. 
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